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ABSTRACT 


After  several  years  of  research,  the  image  theory  has  been 
found  acceptable  to  solve  the  sound  propagation  problem  within 
a  wedge-shaped  fluid  overlying  either  a  slow  or  a  fast  bottom. 
Some  further  observations  where  done  in  this  research.  The 
primary  purpose  was  the  development  and  computer  application 
of  a  new  analysis,  the  doublet  analysis,  based  on  the  image 
theory.  In  this  analysis,  the  features  of  the  image  theory 
are  studied  from  another  point  of  view,  based  on  a  collection 
of  acoustic  doublets. 
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INTRODUCTION 


I  . 

The  propagation  of  sound  in  a  wedge-shaped  ocean  is  a  very 
interesting  problem,  which  is  still  being  investigated.  A 
better  understanding  will  be  significant  for  scientific 
research  and  military  applications. 

Both  the  image  and  parabolic  equation  (PE)  models  have 
been  applied  to  this  problem.  Scientists  who  have  been 
working  with  PE  models  are  listed  in  refs.  1  and  2,  and  those 
with  the  image  theory  in  refs.  3  through  7.  At  the  Naval 
Postgraduate  School  (NPS) ,  recent  theses  on  this  problem,  are: 

1  .  A  PE  model  computer  program  developed  by  Larry  Ernest 
Jaeger  [Ref .2]  . 

2.  Three  other  theses  mostly  about  the  image  theory,  which 
are  listed  in  refs.  5  through  7. 

The  most  recent  thesis  was  done  by  Kim  Jong  Rok,  who,  after 
a  detailed  development  of  the  image  theory,  did  a  comparison 
between  Jaeger's  model  and  image  theory,  producing  results  for 
a  frequency  of  5Hz,a  wedge  angle  of  10°,  a  source  depth  of  50m 
and  various  receiver  depths.  All  analyses  have  shown  that  the 
image  theory  is  acceptable . 

This  research  will  try  to  further  verify  the  reliability 
of  the  image  theory  in  various  cases,  to  identify  its 
behavior,  and  to  develop  an  improved  analysis  based  on  this 
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theory.  More  specifically  this  thesis  will  address  the 
following  tasks: 


1.  Develop  an  image  theory  core  program,  which  will  be  the 
basis  of  any  future  study.  The  associated  computer 
program  will  be  under  the  name  UFTEXT . 

2.  Use  acoustic  doublet  radiation  to  compare  and  analyze 
image  theory  results.  This  comparison  assumes  some 
matching  of  the  water  and  fast  bottom  densities  and  sound 
velocities . 

3.  Verify  the  behavior  of  the  sound  field  at  points  where 
total  internal  reflections  from  the  first  images  start  to 
contribute . 

4.  Extend  a  new  analysis  which  combines  the  images  in  pairs 
of  acoustic  doublets,  and  compare  with  the  original  image 
theory . 
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II. 


DEVELOPMENT 


A.  DEFINITIONS  OF  THE  SYMBOLS 

The  following  symbols  are  used,  (see  Figs.  1  through  9), 
and  all  distances  are  scaled  in  terms  of  the  scaling  distance: 


p 

= 

wedge  angle 

(Eq .  1 ) 

R1 

= 

distance  from  apex  to  source 

(Eq .  9) 

r2 

= 

distance  from  apex  to  the  receiver 

(Eq.  9) 

Y0 

= 

distance  along  the  shore  between  the 

receiver  and  the  source 

(Eq .  9) 

Y 

= 

source  angle  measured  downward 

from  the  surface 

(Eq.7) 

5 

= 

receiver  angle  measured  downward 

from  the  surface 

(Eq .  9) 

Pi 

= 

water  density 

(Eq.  13) 

P  2 

= 

bottom  density 

(Eq  .13) 

Cl 

= 

water  sound  speed 

(Eq  .13) 

c2 

bottom  sound  speed 

(Eq  .13) 

ec 

- 

critical  angle  for  fast  bottom 

(Eq .  1 ) 

e. 

= 

"critical  angle"  for  slow  bottom 

(defined  for  convenience) 

(Eq.4) 

©t 

= 

transmission  angle 

(Eq .13) 

= 

scaling  distance  for  fast  bottom 

(Eq .  1 ) 

scaling  distance  for  slow  bottom 


X 
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(defined  for  convenience) 


ki 

k2 

a/k2 

n 


©„ 


Run 


R 


In 


©„ 


©,. 


0, 


R(©un.) 

R(©lnJ 

Pu 


wavenumber  in  the  wedge  (water) 
wavenumber  in  the  bottom 
loss  term  in  the  bottom 


(Eq.4) 
(Eq.l) 
(Eq .  2) 
(Eq.  17) 


image  number  in  the  description  of 

the  image  theory,  and  pair  of  acoustic 

doublets  number  in  the  description 

of  the  doublet  analysis  (Eq.7,33) 

angle  of  the  nth  image  measured 

from  the  bottom  (Eq.7) 

distance  between  the  receiver  and 

the  nth  upper  image  (Eq.9) 

distance  between  the  receiver  and 

the  nth  lower  image  (Eq.10) 

grazing  angle  of  incidence  on  the  mth 

plane  for  the  nth  image  (Eq.13) 

grazing  angle  of  incidence  on  the  mth 

plane  for  the  nth  upper  image  (Eq.ll) 

grazing  angle  of  incidence  on  the  mth 

plane  for  the  nth  lower  image  (Eq.12) 

reflection  coefficient  corresponding 

to  0unB  (Eq .18) 

reflection  coefficient  corresponding 


to  © 


Inm 


complex  pressure  from  the  upper 
family  of  images 


(Eq.  19) 


(Eq.  18) 
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complex  pressure  from  the  lower 
family  of  images 

distance  between  the  images  of  a 
doublet 

distance  of  the  neutral  doublet 
(the  source  and  its  first  surface 
reflected  image)  from  the  receiver 
receiver-neutral  doublet  axis  angle 
defined  in  page  24 

distance  between  the  receiver  and 
the  nth  upper  doublet 
distance  between  the  receiver  and 
the  nth  lower  doublet 
distance  between  the  receiver  and 
the  nth  doublet,  the  receiver  being 
on  the  surface  (5=0)  for  a  given  R2 
receiver-nth  upper  doublet  axis 
angle 

receiver-nth  lower  doublet  axis 
angle 

receiver-nth  doublet  axis  angle, 
the  receiver  being  on  the  surface 
(6=0)  for  a  given  R2 . 
as  defined  in  eq.45 
as  defined  in  eq.46 
as  defined  in  eq.47 


(Eq  .19) 

(Eq.23) 

(Eq . 31 ) 
(Eq .  31 ) 

(Eq.  33) 

(Eq.34) 

(Eq. 35) 

(Eq.  38) 

(Eq  .39) 

(Eq .  40 ) 
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e 


X'¥ 

z 

Xn+/¥n+ 

Xn~  f  ¥n~ 

P„+ 

Pn“ 


as  defined  in  eq.48 

as  defined  in  eq.49 

represent  the  different  reflection 

coefficients  of  the  images  of  an 

acoustic  doublet  (Eq.25) 

combined  reflection  coefficient 

for  a  doublet  (Eq.26) 

reflection  coefficients  for  the  nth 

upper  doublet  (Eq.52) 

reflection  coefficients  for  the  nth 

lower  doublet  (Eq.52) 

complex  pressure  for  the  nth  upper 

doublet  (Eq.43) 

complex  pressure  for  the  nth  lower 


doublet  (Eq.44) 

Pn  =  complex  pressure  for  the  nth  pair  of 

acoustic  doublets  (Eq.50) 

It  is  worth  noting  that  the  source  angle  y  and  the 
receiver  angle  8  are  measured  downward  from  the  surface, 
although  in  previous  theses  they  were  measured  upward  from  the 
bottom.  It  was  found  more  convenient  to  do  so  for  the 
purposes  of  this  research. 


B .  IMAGE  THEORY 

The  theoretical  background  of  the  image  theory  has  been 
presented  by  Coppens,  Sanders,  Ioannou  and  F.awamura  [Ref. 3]. 


This  section  will  be  a  brief  review  of  the  ideas  and  basic 
equations  of  the  image  theory. 

The  assumptions  are: 


1 .  The  sound  velocities  and  densities  are  constant  both  in 
the  wedge  and  in  the  penetrable  bottom. 

2.  The  water-air  boundary  is  a  pressure  release  surface. 

3.  The  slope  is  constant. 

4.  All  boundaries  are  planar. 


All  distances  are  scaled,  according  to  R.cal<Kl  =  Rr.«i/X.  For 
a  fast  bottom,  the  scaling  distance  Xc  is  defined  as  the 
distance  where  the  lowest  mode  attains  cut  off.  For  a  slow 
bottom,  the  lowest  mode  does  not  attain  cut  off,  but  a  scaling 
distance  X.  is  defined  for  convenience  in  calculations. 

The  definitions  are 
Fast  bottom 


klX°  2  sin0ctan|i 
kiXc  2  tanOctan|l 


cos  9  = 

C2 


Slow  bottom 

]e  X  -  _ — _ 

1  '  i  tan9#tan[J 


(1) 

(2) 

(3) 

(4) 
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where  (3  is  the  wedge  angle,  subscript  1.  refers  to  the  fluid  in 
the  wedge  (water) ,  and  subscript  2  refers  to  the  bottom. 

The  number  of  images,  including  the  source,  which 
contribute  to  the  sound  field,  is  given  by  n  =  INT(360°/P), 
(Fig.l) .  The  angles  between  each  image  and  the  bottom  (Figs. 
1  through  3)  are  given  by 

0n  =  nP-y,  for  n  odd  (7) 

@n  =  (n-l)p^y,  for  n  even  (8) 

where  y  is  the  angle  of  the  source  measured  from  the  surface . 

The  distance  between  each  image  and  the  receiver  (Figs .  2 
and  4) ,  is  given  by 


Ran  =  V*i2 -2*^003  (0„-P  +  5)  -F02 

for  the  upper  group  of  images,  and 

Rln  =  ^^22-2^2cos  (0„  +  p-5)  +3To 


(9) 


(10) 


for  the  lower  group  of  images,  where  Rx  is  the  source-apex 
distance,  R2  is  the  receiver-apex  distance,  5  is  the  angle  of 
the  receiver  from  the  surface,  and  Y0  is  the  horizontal 
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Figure  1.  Image  Structure  for  a  Wedge-shaped  Duct 
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source  plane 


receiver  plane 


bottom  plane 


Figure  2.  Illustration  of  the  Angles  ©n,  y,  6,  and  the  Angle  between 
nth  Image  Plane  and  the  Receiver  Plane. 
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Figure  3.  The  Angle  of  Incidence  0un„  of  the  nth  Upper  Image  on  the  mth 
Plane . 


Figure  4.  The  Distance  R^, 


source-receiver  distance  along  the  shore.  All  distances  are 
scaled . 

The  incident  angle  of  the  nth  image  to  the  mth  plane  (Figs . 
1  and  3)  ,  is  given  by 


.  Kxsin  (@n-2mP)  +R2  sin  [  (2m-l)  p  +  S]  (11) 

sin  (v?  )  - — - r  m-1 , 2 ,  .  .  . 


for  the  upper  group  of  images,  and 


•  i  -  -RiSin(0n-2mP) +R2sin[  (2m+l)p-6]  (12) 

sm  (“iIU# )  - — - ,  m-0 ,  1 ,  .  .  . 


'In 


for  the  lower  group  of  images . 

The  reflection  coefficients  are  assumed  to  be  the  plane 
wave  Rayleigh  reflection  coefficients,  a  good  approximation 
for  the  problem.  The  equations  for  them  are: 

Without  absorption 

sin  (0^)  -Zisin  (0t) 

«(©-)  =  - £ -  (13) 

^isin(0Jun)  + — _  sin  (0t) 

Pi  C2 


sin©t 


N 


1  _  ( — ) 2  cos2  <©„) 
C1 


(14) 


where  0^,  is  the  incident  grazing  angle  of  the  nth  image  on  the 
mth  plane  (Fig.  3)  and  ©t  is  the  corresponding  transmission 


grazing  angle. 

With  absorption 
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(15) 


a  =  {c2/ c2)  2-cos20m  (16) 

Jb  =  2  {cjc2)2a/k2  (17) 

where  a  is  the  absorption  coefficient. 

The  total  sound  field  pressure  contribution  from  the  upper 
images  at  each  point  in  the  wedge,  is  given  by 

Pu  =  £-Lexp(-j ^Run)  (-l)'"/2)nR(0UM)  (18) 

riR(0unin)  =  1,  for  n  -  1,2  and  M  =  (n-l)/2, 
and  from  the  lower  images, 

P 1  =  ri-exp(-#1Rln)  (-l)'"/2,nR(0Im)  (19) 

-=0 

The  total  pressure  will  be,  P  =  Pu+Px. 

C .  ACOUSTIC  DOUBLET  RADIATION  WITHOUT  REFLECTION  COEFFICIENTS 

An  acoustic  doublet  consists  of  two  simple  sources  of 
equal  strengths,  separated  by  a  distance  d,  vibrating  with  the 
same  frequency  but  180°  out  of  phase  with  each  other  [Ref. 8] . 

Assume  as  a  reference  that  the  source  closer  to  the 
receiver  is  the  negative  one  and  the  other  one  is  the  positive 


14 


(Fig.  5)  .  The  pressure  at  point  (r,  0)  due  to  the  positive 
source  is 


Pl  T^77e 


and  that  due  to  the  negative  source  is 


J[att  -k  (r-Ai-j ) ) 


(21) 


where  r  is  the  distance  from  the  field  point  to  the  midpoint 
between  the  sources,  Arx  and  Ar2  are  as  shown,  and  the  opposite 
signs  account  for  the  phase  difference  of  the  sources.  The 
acoustic  doublet  pressure  at  (r,0)  is  then 


(22) 


The  process  of  replacing  Arx  and  Ar2  by  r  and  0  is 
straightforward  but  tedious,  and  the  resulting  expression  for 
p(r,0,t)  is  too  complicated  for  efficient  analysis.  However, 
in  the  most  frequently  encountered  cases,  as  it  happens  in  the 
wedge-shaped  ocean  problem,  the  observation  of  the  pressure  is 
made  at  distances  great  compared  to  the  separation  of  the 
sources.  Therefore,  it  will  be  useful  to  derive  the  form  of 
eq.22  suitable  in  the  limit  r»d.  This  is  referred  to  as  the 
far-field  approximation.  If  r>>d,  Arj«Ar2=»  (d/2)  sin0  and 
Ar1/r<<l,  so  that  eq.22  becomes 

p(r,  0,t)  =  -  j-^sin  ( -i)cdsin0)  ej{a>t  kr)  (23) 
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figure  5.  Geometry  Used  in  Deriving  the  Radiation  Characteristics  of  an 
Acoustic  Doublet,  according  to  Eq.23. 
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If  the  signs  of  the  sources  are  reversed  (Fig. 6),  then 


eq.23  changes  sign  too,  and  becomes 

p(r,  8,  t)  =  j-H^sin  (.IjcdsinQ)  eji<ot'kr)  (24) 


D.  ACOUSTIC  DOUBLET  RADIATION  ENCOUNTERING  REFLECTION 
COEFFICIENTS 

Consider  again  the  situation  of  Fig. 5,  which  was  analyzed 
in  the  previous  section.  Suppose  now  that  the  free  acoustic 
fields  of  the  positive  and  negative  sources  encounter 
reflection  coefficients  x  and  V  respectively.  Those 
reflection  coefficients  are  assigned  to  correspond  to  grazing 
angles  0X,  0V  as  shown  in  Fig. 7.  Then  eq.22  becomes 


P  = 


A 

r 


e-J*Ar. 

1  +Ar1/r 


~^KF2 Ti 


)  eJ(o>t-*r) 


(25) 


Combination  of  eq.23  and  eq.25  gives 


A 

r 


(X 


'  j~2 


1  **d«±n0 


^  -\\t  w  )  = 

1  +Ar1/ r  T  1  -Ar2/  r 


-j2—z  sin  ( .ifcdsin  0) 
r  2 


(26) 


for  the  far  field  approximation  (r»d)  ,  where 
Ar1=Ar2=  (1/2)  dsin0  and  (Ari/r)<<l,  (Ar2/r)  «1 .  The  coefficient 
.z  is  the  appropriate  reflection  coefficient  for  the  acoustic 
doublet  field.  After  some  mathematical  manipulation  of  eq.26, 


z  =  <X+V) -j(V-X)cot  (_^/cdsin0)  ] 


(27) 
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Figure  6.  Geometry  Used  in  Deriving  the  Radiation  Characteristics  of  an 
Acoustic  Doublet  according  to  Eq.24. 
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Figure  7.  Geometry  Used  in  Deriving  the  Radiation  Characteristics  of  an 
Acoustic  Doublet,  Including  Reflection  Coefficients,  according  to  Eqs . 
25  and  26. 
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Suppose  now  that  the  polarity  of  the  sources  changes 
(Fig. 8),  and  consider  eq.24.  In  that  case,  the  signs  in  both 
sides  of  eq.26  change,  and  the  solution  of  eq.27  for  the 
doublet  reflection  coefficient  z  remains  the  same. 

Surprisingly,  this  reflection  coefficient  didn't  turn  out 
to  be  too  complicated  to  apply  in  the  doublet  analysis  of  the 
image  theory,  because  it  is  expressed  as  a  function  of  the  sum 
and  the  difference  of  the  reflection  coefficients  of  the 
individual  images  in  the  doublet . 

In  this  research,  the  choice  for  the  difference  was  the 
reflection  coefficient  of  the  image  closer  to  the  receiver 
minus  the  reflection  coefficient  of  the  image  further  from  the 
receiver.  That  can  be  seen  from  eq.27  in  connection  with  Figs 
7  and  8 . 

Combination  of  eq.27  with  eqs .  23,  24  and  26,  gives  for 
the  two  different  polarities  the  total  acoustic  doublet  fields 

p  =  [  <Y~X)  cos  ( -i-kdsinO)  +  j  (X+Y)  sin  ( -^JcdsinO)  ]  e  ) 

p*  =  [  (Y~X)  cos  ( -i/cdsin9)  +  j  (X+Y)  sin  ( -ifcdsinO)  ]  ej{at  kr)  ^  ^ 

with  A  has  been  set  equal  to  one . 

These  two  equations  were  successfully  applied  in  the 
doublet  analysis  of  the  image  theory  as  will  be  discussed 
later . 
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Figure  8.  Geometry  Used  in  Deriving  the  Radiation  Characteristics  of  an 
Acoustic  Doublet,  Including  Reflection  Coefficients,  according  to  Eq.29. 
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E.  DOWNSLOPE /UP SLOPE  DOUBLET  ANALYSIS  OF  THE  IMAGE  THEORY 


1 .  General  Description 

The  objective  of  this  analysis  is  to  combine  the 
images  in  pairs  of  acoustic  doublets  and  thus  simplify  the 
image  theory  [Ref. 9]. 

The  procedure  starts  (Fig. 9)  by  first  considering  the 
1**  and  2nd  upper  images,  the  source  and  its  first  surface 
reflected  image,  as  a  neutral  doublet.  It  is  obvious  that 
this  neutral  doublet  (n=0)  does  not  encounter  any  reflection 
coefficients . 

The  3rd  and  4th  upper  images  (Fig.  9)  constitute  the  I*1 
upper  doublet,  and  the  l'*1  and  2nd  lower  images  constitute  the 
lBt  lower  doublet .  Both  doublets  constitute  the  l't  pair  of 
acoustic  doublets  (n=l) .  This  procedure  continues  until  all 
the  images  are  organized  in  pairs  of  acoustic  doublets. 

The  total  number  of  pairs  of  acoustic  doublets  is 


N  =  Int 


[Int  (2^1)  -2] 
- 3 - 


(30) 


It  is  easily  seen  that  in  the  worst  case  three  images 
might  be  missed.  In  that  case,  those  will  be  images  at  the 
end  opposite  to  the  source.  Their  contribution  will  be 
negligible,  because  many  reflection  coefficients  will 
contribute  to  them,  and  among  those  reflection  coefficients 


t 
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figure  9.  Geometry  Used  in  Deriving  the  Doublet  Analysis. 
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will  be  some  with  very  small  values.  Hence,  those  images  will 
be  ignored. 

2 .  Neutral  Doublet  Acoustic  Field 

The  neutral  doublet,  which  consists  of  the  source  and 
its  first  surface  reflected  image  (Fig. 9),  does  not  encounter 
any  reflection  coefficients. 

Application  of  eq.24  gives  the  pressure  from  the 
neutral  doublet, 

P0  *  +j_!sin  (i.^dsinOoJe^'*^’  (31) 

ro  *- 

~  -  =  k,R, Y 

sinO0  =  — smo  =  — o 


1  R.  R,  _ 

^./CjdsmOpi  k1 - yo 

—  n 


If  we  define 


R  R7 
r  -  k  12 

1  0  A.x. 


then  the  pressure  has  the  simple  form 


P0*  +  j_fL  sin  <roy6)  eJ,(a,t'*‘r°) 


(32) 


where  P0  is  the  neutral  doublet  complex  pressure,  r0  is  its 
distance  from  the  receiver,  and  CT0  is  the  receiver-doublet- 


doublet  axis  angle.  The  positive  sign  in  front  of  the  right 
hand  side  of  eq.31  and  eq.32  counts  for  the  positive  reference 
of  the  source  (Fig. 9) . 

3 .  Acoustic  Field  from  the  Pairs  of  Acoustic  Doublets 
Ignoring  the  Bottom  Reflection  Coefficients 

The  first  pair  of  acoustic  doublets  is  shown  in  Fig. 9. 
The  distances  from  the  receiver  of  the  individual  upper  and 
lower  doublets  at  2n|3  are 


rn*  =  /r12+R22-2R1J?2cos  (2np  +  5) 
for  the  upper  doublet  and 


(33) 


r 


n 


*R|-2R1R2cos  (2n(3-5) 


(34) 


for  the  lower  doublet. 

When  the  receiver  is  at  the  surface,  then  the  above 
distances  are  equal 


r 


nO 


12+R22-2JR1R2cos  (2 np) 


(35) 


Eqs .  33 
Taylor's  series 


and  34  to  the  2nd  order  approximation  of  the 
expansion  around  8=0,  are 


_  ■  _  ^rn 

"  r”°  °-35‘ 


6*0 


-is2  d 

2 


rr 

dV 


6*0 


drn‘ 

"38" 


6«0 


— i —  sin  (2 n(3) 
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dr. 


^}—l  sin  (2 n{3) 


5=0  n0 


d2r 


d82 


5=o 


d2rn 

"W 


RiR 

6=0  ^nO 


[cos  (2n(5)  -  _ sin2  (2np)  ] 


rn0 


-r„*  =  rn0+J^i8  [sin  (2n|5)  +  -i.8cos  (2n(5)  -_i &i!l^isin2  (2np)  ] 

rn0  Z  Z  r„o 


1  X  ^1^2  _  ^ 

r20 


(36) 


rn  =  rn0^5...i  5  [  -  sin  (2nf3)  +  _i.8cos  (2n|3)  -.1.8 — L_isin2  (2nP) 

2  2  r~  ~ 

n 


’  nO 


•  nO 


Receiver-doublet-doublet  axis  angles  are 


s^r  o„*)  =  isin(2np  +  5) 


for  the  upper  doublet  and 


sin(an  )  --  — Lsin(2n|3-8) 
r  ■ 


(38) 


(39) 


for  the  lower  doublet . 

When  the  receiver  is  at  the  surface  this  angle  becomes 
the  same  for  both  doublets. 


sin(GnC))  =  — isin(2n|3) 

rnO 


(40) 


Eqs .  38  and  39  to  the  l-t  order  approximation  of  the 

Taylor's  series  expansion  around  8=0  are 
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sin  (<J„)  =  sin  (CTn0)  +8 


dsin. 


SF 


8=o 


dsin  (o') 


dd 


R  R  R 

_ 1  [cos  (2n(5)  -  i  2  sin2  (2n|J)  ] 


6  =  0  An0 


-^rO 


=  — icos  (2nB) 
r_„ 


dsin (cr  ) 


d& 


R  R  D 

— L  [cos  (2 np)  - — i_lsin2  (2n(5)  ] 


8=o 


’  n0 


2 

*^*n0 


cos  (2np) 


sin(On‘)  =  — L  [sin  (2n{5)  +6cos  (2n^)  ] 


(41) 


sin  (cvn  )  =  — i  [ sin  (2n(5)  -8cos  (2nf3)  ] 


(42) 


The  complex  pressure  field  for  the  nth  upper  doublet  is 


pn'  =  ("1)n  j-~  sin  [-i^dsin  <C„‘)  ]  e 

IT_  ^ 


(43) 


and  for  the  n  lower  doublet 


pn  =  -(-l)n  j—  sin  [—k^dsin  (On~)  ]  e 
r  '  2 


jtcot-k,  r.) 


(44) 


The  approximation  (l/2)k1d=R1y  and  the  definition 


r  =  k. 


r,r2 


(45) 


'  n0 


give 


=  rnysin(2nP)  =  -i/qdsin  (an0) 


(46) 
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<p„  =  rn6sin(2np)  =  kxb 


dr. 


~d5 


8=o 


(47) 


dsina. 


£„  =  rny6cos(2nP)  =  (^k^d)  8 — ^ 


8=o 


(48) 


„  1  §£  (p  2  i 

E  =  jL  ( _ 71  -  Yn  ?  =  _£*. 

“n  ^  L-  —  '  1 


d2r 


*ir„c 


dS5 


6=o 


(49) 


The  above  approximations  and  the  definitions  in  eqs . 
45  through  49,  give  the  pressure  field  for  the  nth  pair  of 
acoustic  doublets  (the  sum  of  the  pressure  fields  of  the 
corresponding  upper  and  lower  doublets) 


pn  =  pn+pn  =  ( - 1 )  "  jJL  e j{wt [  e  "*■  sin  (d„  *  E  „)  -  e *■  s  in  ( $n - E  „)  ] 

r"°  (50) 

More  specifically,  the  following  approximations  where 
applied  in  eq.50: 

1.  rn0  used  in  amplitude  (2/rn0),  instead  of  rn+  and  rn_ . 

2.  rn  correct  to  the  2nd  order  approximation  (2;n)  in  phase. 

3.  <Tn  correct  to  the  l*t  order  approximation  (En)  in 
amplitude . 

The  total  field  is  the  sum  over  all  pairs  of  acoustic 
doublets,  plus  the  field  of  the  neutral  doublet: 

N 

p  =  p0-£  lPn‘+pB-]  <51> 
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4 .  Acoustic  Field  from  the  Pairs  of  Acoustic  Doublets 
Encountering  the  Bottom  Reflection  Coefficients. 

The  reflection  coefficient  of  eq.27  is  applied  to 
eq.50.  Hence,  according  to  eqs .  28  and  29,  the  pressure  field 
including  the  bottom  reflection  coefficients  is 


P„  =  (-1)”  _LeJ(,Dt-**r-' u(e -*•[  (cos(Vg  )  (Y^-Xn4) +7  (sin  (0n  +  en) ) 


(52) 

(X/+V„‘)  ]  -«*•[  (cos  (0„-£n)  )  (Vn  -X„‘)  +  J  (sin  (fln-£„)  )  (x„‘  +Vn')  1* 

where  \|/n+,  \j/n~  are  the  reflection  coefficients  of  the  upper  and 
lower  doublet  images  closer  to  the  receiver,  and  Xn  r  Xn  are 
those  further  from  the  receiver.  The  corresponding  grazing 
angles  0X+,  0X",  0V+,  0V~  are  illustrated  in  Fig.  9. 

The  total  field  is  given  by  eq.51,  where  Pn=Pn++Pn~  is 
given  by  eq.52. 
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Ill . 


COMPUTER  PROGRAMS 


The  following  computer  programs,  were  created  and  used 
during  this  research.  They  were  written  in  the  FORTRAN  77 
computer  language,  on  the  IBM-3033  main  frame  NPS  computer 
system.  Some  of  the  programs  use  plotting  routines,  using  the 
DISSPLA  graphics  package,  while  URTEXT  and  DOUBLT1C  do  not. 
Copies  of  all  programs  are  available  in  professors  Coppens 
main  frame  account . 

A .  PROGRAM  URTEXT 

In  the  beginning,  the  program  URTEXT  was  created  to 
provide  a  well  documented,  immutable  program  of  the  wedge- 
shaped  ocean  image  theory.  The  follow-on  programs  are 
modifications  of  URTEXT. 

In  URTEXT,  the  receiver,  at  distance  R2  from  the  apex, 
moves  from  the  surface  to  the  bottom  of  the  wedge,  on  a 
circular  arc  (Fig. 10) .  This  program  provides  two  distinct 
results  of  the  image  theory: 

1.  Pressure  amplitude  vs  receiver  angle. 

2.  Phase  vs  receiver  angle. 

The  receiver  angle  is  measured  downward  from  the  surface. 
This  program  accepts  input  data  and  also  gives  the  output  data 
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with  respect  to  the  scaling  distance,  independent  of  the 
frequency.  Appendix  "B"  is  a  copy  of  URTEXT. 

B .  PROGRAM  URTEXT 1 

In  this  program  the  receiver  moves  horizontally.  Hence, 
its  depth  from  the  surface  is  kept  constant  (Fig. 11).  This 
program  provides  two  distinct  results  for  comparison: 

1.  Transmission  loss  (TL)  vs  source-receiver  horizontal 
range,  using  the  complete  image  theory. 

2.  Transmission  loss  (TL)  vs  source-receiver  horizontal 
range,  using  only  the  neutral  doublet  acoustic  radiation. 

The  inputs  and  outputs  are  frequency  dependent  and 
distances  are  expressed  in  meters  or  kilometers. 

C .  PROGRAM  URTEXT 1C 

In  this  program,  the  receiver  moves  as  does  in  URTEXT  and 
it  provides  the  same  results.  The  only  difference  is  the 
addition  of  a  plotting  routine,  to  present  the  data. 

D .  PROGRAM  D0UBLT1C 

In  this  program,  the  receiver  moves  the  same  as  in  URTEXT, 
and  it  also  provides  the  same  results .  The  calculations  are 
done  according  to  the  doublet  analysis  presented  previously, 
without  including  any  absorption  in  the  bottom  (ot/k2=0)  or  in 
the  water.  The  reflection  coefficient  used,  was  that  of 
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Figure  11.  Positions  of  the  Receiver  Considered  in  Program  URTEXT1 
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equation  13.  It  works  for  a  slow  bottom  only.  It  was  used 
with  their  real  reflection  coefficients,  this  program  was  used 
to  simplify  the  initial  verification  of  the  doublet  analysis. 
Appendix  "C"  is  a  copy  of  this  program. 

E .  PROGRAM  CALI 

This  program  calculates  the  cosine  of  the  incident  grazing 
angle  of  the  nth  image  on  the  mth  plane,  as  a  function  of 
range.  It  was  used  to  determine  at  what  range  from  the  source 
the  total  reflection  contribution  from  the  first  images 
starts,  as  a  function  of  the  critical  angle,  defined  as 

COsO^Ci/Cj  . 

F.  PROGRAM  PLOTIFD 

This  is  a  plotting  program  developed  by  Jaeger  to  plot  his 
PE  model's  data.  It  was  modified  to  overplot  URTEXT1  data  for 
comparison . 


IV. 


VERIFICATION  OF  THE  PROBLEMS 


A.  COMPARISON  BETWEEN  THE  IMAGE  THEORY  AND  JAEGER'S  PE  MODEL 

WITH  RESPECT  TO  THE  ACOUSTIC  DOUBLET  RADIATION. 

This  was  done  by  using  the  URTEXT1,  JAEGER  and  PLOTIFD 
programs . 

1 .  Cases  Tested 

The  cases  tested  (Figs.  12  to  22)  were: 

1.  Wedge  angles  p  =  10°  and  p  =  5° . 

2 .  Frequency  f  =  50  Hz . 

3.  Source  depth  =  50  m. 

4.  Source  angles  y  =  2.4952°,  2.5048°,  5.0384°. 

5.  Receiver  depth  =  75  m. 

6.  D1  =  px/p2  =  0.9,  0.94,  0.9999. 

7.  CC  =  c1/c2  =  0.9,  0.94,  0.9999. 

8.  Source-receiver  horizontal  ranges  from  0-1  km  and  from  0- 
10  km,  downslope. 

In  these  figures,  the  dashed  line  represents  the  image 
theory  data,  the  dotted  line  represents  the  doublet  (source 
and  first  surface  reflected  image)  data,  and  the  solid  line 
represents  the  Jaeger's  PE  model  data. 

2 .  Specific  Observations 

1.  In  close  distances  from  0-200  m  (Figs.  12  and  13),  the 
image  theory  is  almost  in  agreement  with  the  acoustic 
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doublet  radiation  curve,  while  the  PE  model  does  not 
agree  very  well. 


2.  At  larger  distances,  it  is  hard  to  say  if  the  image 
theory  or  PE  model  is  better,  because  the  contribution  of 
the  images  becomes  significant.  More  specifically: 

a.  The  two  TL  curves  seem  to  agree  with  each  other, 
being  far  away  from  that  of  the  acoustic  doublet. 

b.  In  almost  identical  matching,  D1=0.9999,  CC=0.9999 

and  at  {5=10°  {Fig. 21),  all  three  curves  were 
perfectly  coincident.  The  same  case  at  f$=5° 
(Fig. 22),  gave  a  difference  between  the  image  theory 
and  the  acoustic  doublet,  starting  from  about  3  km, 
and  reaching  about  2.5  db  at  10.0  km,  but  gave  a  good 
coincidence  up  to  3  km.  This  might  have  happened  due 
to  the  more  significant  contribution  of  more  images 
at  larger  distances . 

3.  For  p=10°  and  D1=0.9,  CC=0 . 9  or  closer  to  unity,  for 
short  ranges  up  to  100  m,  the  effects  of  the  bottom  can 
be  ignored  with  very  small  error. 

B.  CONTRIBUTION  OF  THE  IMAGES 

By  using  the  program  CALI,  some  calculations  were  made  of 
the  grazing  angle  of  the  nth  image  with  the  mth  plane  versus 
range.  When  this  angle  is  equal  to  the  0C,  this  is  the  range 
where  total  internal  reflection  of  the  image  starts.  This 
test  was  made  for  the  l'*  and  2nd  lower  images  and  the  3rd  upper 
image,  because  these  contribute  most  significantly.  The 
following  observations,  were  done  for  the  cases  considered 
previously : 


1.  At  (5=10°,  CC=0.9,  the  l"1  lower  image  contributes  total 
reflection  at  236  m,  the  2nd  lower  image  at  533  m,  and  the 
3rd  upper  image  at  764  m.  Almost  at  those  locations  three 
characteristic  bumps  appear  on  the  image  theory  TL  curve 
for  that  case  (Fig. 13) .  On  the  PE  model  curve  the  first 
bump  appears  correctly,  the  second  did  not  appear,  and 
the  third  appeares  with  a  difference  of  33  m  (Fig. 13). 
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2.  The  same  case  was  tested  at  (5=5°,  CC=0.9,  and  the 

distances  obtained  where  187  m, 427  m,  581  m.  The  image 
theory  TL  curve's  bumps  correctly  appeared  almost  at 
those  distances  (Fig. 15). 


All  tests  done  above  show  that  the  image  theory  behaves  as 
expected.  Hence,  it  is  reconfirmed  that  it  is  deemed  worthy 
to  be  applied  in  the  wedge-shaped  ocean  problem. 

C.  COMPARISON  BETWEEN  THE  IMAGE  THEORY  AND  THE  DOUBLET 
ANALYSIS 

This  comparison  was  done  using  the  results  of  URTEXT1C  and 
DOUBLT1C .  Some  of  the  results  obtained  can  be  seen  in  Tables 
1  to  5,  in  Appendix  "A".  Those  tables  clearly  show  the 
excellent  agreement  between  the  two  programs.  Many  other 
cases  were  tested  to  verify  those  results,  but  are  not 
included  in  this  thesis.  All  results  that  were  tested  were  in 
agreement.  More  specifically: 


1 .  In  all  cases  tested,  there  was  an  excellent  agreement  in 
amplitude,  as  it  can  be  seen  from  some  data  collected  in 
Tables  1  to  5. 

2.  Although  very  good  agreement  (3rd  decimal)  is  still  valid 
in  amplitude  when  the  source  gets  close  to  the  bottom, 
in  that  case  appears  a  phase  difference  not  larger  than 
30°  in  most  cases  tested.  An  example  can  be  seen  in  the 
data  collected  in  Table  1 .  The  reduced  accuracy  of  the 
doublet  approximation  (Eqs  .  23  and  24),  due  to  the  larger 
d,  is  not  responsible,  because  the  phase  difference  does 
not  disappear  even  when  RJ-R1>>d.  It  was  seen  that  with 
higher  order  approximations  of  the  Taylor' s  series 
expansion  of  rn,  an  in  both  phase  and  amplitude  of  eq.52, 
that  phase  difference  reduces.  Probably,  it  is  not  worth 
complicating  the  analysis  by  including  these  higher 
orders,  because  the  accuracy  as  it  exists  will  be 
acceptable  in  most  practical  cases .  Also  note  that  only 
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the  amplitude  is  necessary  for  the  calculation  of  the 
transmission  loss,  and  that  is  obtained  with  excellent 
accuracy . 

Because  of  the  good  agreement  described  previously,  the 
doublet  analysis  can  be  used  as  a  substitute  of  the  full  image 
theory . 

During  this  research,  efforts  were  made  to  simplify  the 
doublet  analysis.  It  was  found  that  the  approximations 
described  in  sect.  E  of  the  development,  can  not  be  simplified 
to  give  correct  results  everywhere  in  the  wedge. 

Consider  the  grazing  angles  formed  between  rn+  and  bottom 
planes  for  the  upper  doublets,  and  between  rn"  and  bottom 
planes  for  the  lower  doublets  (Fig. 9).  A  further 
simplification  of  the  doublet  reflection  coefficients  could  be 
accomplished  by  considering  that  the  reflection  coefficients 
correspond  to  these  grazing  angles.  This  method  gives  correct 
results  only  under  the  follow  conditions: 

1.  Large  wedge  angles  (P>5°),  which  usually  are  not  met  in 
practice . 

2 .  The  receiver  shallower  than  the  surface-bottom  middle 
depth . 

3.  Rx/R2<<1  cannot  hold  for  the  downslope  case,  nor  R1/R2>>1, 
for  the  upslope  case. 

It  was  not  possible  to  determine  an  exact  threshold  to 
apply  such  an  approximation,  because  the  total  field  is 
calculated  encountering  multiple  factors.  For  that  reason, 
the  reflection  coefficient  as  given  in  eq.27  was  used.  Some 
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explanations  of  the  necessity  to  use  eq.27  for  the  reflection 
coefficient,  can  be  derived  by  examining  eq.28.  This  equation 
is  reproduced  here  without  the  phase  factor  e3(®t"kr)  : 

p~  =  --j  [  (ty -X)  cos  ( -i/cdsin©)  +j  (x  +  i|f )  sin  {  —  kdsinQ)  ] 

X+V  is  always  greater  than  \y~x,  for  real  reflection 
coefficients.  The  smaller  the  difference  between  the 
reflection  coefficients  x>¥  of  the  images  in  the  doublet,  the 
more  dominant  the  imaginary  part  of  the  above  equation  will 
be,  with  respect  to  X'V- 

Consider  now,  small  wedge  angles  ((5<5°),  along  with  R1«R2 
for  downslope  and  Ri»R2  for  upslope.  Due  to  the  small  wedge 
angle,  the  lower  images  contribute  significantly.  Then,  for 
these  lower  images,  cos  [  (1/2)  kdsinO]— >1  and  sin  [  ( 1/2)  kdsin9]-*0 , 
because  of  the  small  value  of  the  receiver-doublet  axis  angle, 
under  the  above  conditions.  Therefore,  in  that  case,  the  real 
part  will  be  more  dominant. 

It  becomes  clear  from  the  above  discussion,  that  in  many 
cases,  neither  the  real  or  the  imaginary  parts  of  eq.28  can  be 
ignored.  Thus,  eq.23  is  the  correct  equation  to  apply  for  the 
wedge . 
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Figura  13.  Downslope  Transmission  Loss  (dB)  vs  Horizontal  Source- 
Receiver  Distance,  to  1km,  for  JJ-10®,  7-5.0384*,  Pi/p2-0.9,  C,/C2-0.9. 
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HORIZONTAL  RANGE  (KM 


Figure  15.  Downslope  Transmission  Loss  (dB)  vs  Horizontal  Source- 
Receiver  Distance,  to  1km,  for  £-5®,  7-2.5048®,  Pi/pj-0.9,  Ci/C;-0.9. 
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Figure  17.  Downslope  Transmission  Loss  (dB)  vs  Horirontal  Source- 
Receiver  Distance,  to  1km,  for  P“10*,  7*5. 0384*,  pi/Ps^O . 94,  Cj/Cj-0.94. 
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rigr uta  20.  Downslope  Transmission  Loss  (dB)  vs  Horizontal  Source 
Receiver  Distance,  to  10km,  for  |J“5®,  y— 2.5048*,  Pi/Pj“0.94,  C,/C2""0.94. 
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figuxt  21.  Downslope  Transmission  Loss  (dB)  vs  Horizontal  Source 
Receiver  Distance,  to  10km,  for  P-10’,  y-5.0384*,  Pj/pj-0 . 9999, 

C,/C,-0. 9999. 
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V. 


CONCLUSIONS  AND  RECOMMENDATIONS 


Image  theory  is  commonly  accepted  as  a  solution  to  the 
wedge-shaped  ocean  acoustic  propagation  problem,  under  the 
assumptions  defined  in  sect.  B  of  the  development.  Doublet 
analysis  is  a  successful  substitute  for  the  image  theory. 
This  is  important  for  two  reasons: 


1.  It  is  a  more  compact  solution. 

2.  The  corresponding  computer  program  execution  time  is 
faster  than  that  of  the  image  theory,  because  it  requires 
fewer  loops  to  calculate.  Some  tests  which  were  done 
have  shown  that  tD=0.55tj,  where  tD  is  the  doublet 
analysis  computer  program  execution  time  and  tj  is  the 
image  theory  analysis  computer  program  execution  time. 


It  is  recommended  that  the  doublet  analysis  be  extended 
for  the  cross  slope  case,  and  that  the  corresponding  computer 
program  be  extended  for  a  fast  bottom  and  to  include 
absorption  in  the  bottom. 

It  also  appears  possible  an  additional  simplification,  can 
be  made  according  to  the  dipole  approximation  [Ref. 8],  where 
the  doublet  acoustic  pressure  without  including  any  reflection 
coefficients 

p(r,0,t)  =  - j'^^sin  (-i./c,dsin9)  e^(<Dt  k,r> 
r  2 

further  reduces  to 
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p(r,e,t)  -  ~j— i^sinee:i(a>t-*’r, 

r 

when  k1d=2k1R1y«l .  Using  equations  1  and  4,  this  can  also  be 
expressed  with  respect  to  the  scaling  distance,  as 


for  a  fast  bottom,  and 


for  a  slow  bottom.  It  is  obvious  from  the  last  two  equations 
that  the  dipole  approximation  can  be  applied  whenever  Rj«X . 
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APPENDIX  "A 


NUMERICAL  RESULTS  OF  COMPARISON  BETWEEN  THE  IMAGE  THEORY  AND 
THE  DOUBLET  ANALYSIS 

TABLE  1 

Data  for  B=7  °  ,  V=4  °  ,  R,=4 ,  R,=25 ,  O,/p,=0 . 7 ,  c,  =  1500m/s ,  c,=1430m/s 


PRESSURE 

AMPLITUDE 

PHASE 

ANGLE 

RECEIVER 

ANGLE 

IMAGE 

DOUBLET 

IMAGE 

DOUBLET 

THEORY 

ANALYSIS 

THEORY 

ANALYSIS 

o 

r~ 

o 

o 

0.005314 

0.005282 

58.0° 

85.4° 

o 

o 

0 

0.055618 

0.055174 

56.9° 

84 .4° 

1.40° 

0.124116 

0.123650 

54.4° 

81.8° 

2 . 10° 

0.212473 

0.212157 

51.3° 

78.6° 

2 . 80° 

0.311260 

0.311706 

48.7° 

75.8° 

3.50° 

0.394972 

0.396090 

46.8° 

73.7° 

4.20° 

0.434343 

0.435430 

45.4° 

72.3° 

4 . 90° 

0.412928 

0.413046 

44.1° 

70.9° 

5.60° 

0.331062 

0.329867 

41 . 6° 

68.3° 

6.30° 

0.200917 

0 . 198261 

33.5' 

59.6° 

7 . 00° 

0.069522 

0.072188 

-25.4° 

-4.0° 
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TABLE  2 


Data  for  B=6°  ,  7=2°  ,  R,=0 . 8 ,  R,=5,  D,/o,=0 . 8,  c-,=1500m/s,  c,=1450m/s 


RECEIVER 

ANGLE 

PRESSURE  AMPLITUDE 

PHASE  ANGLE 

IMAGE 

THEORY 

DOUBLET 

ANALYSIS 

IMAGE 

THEORY 

DOUBLET 

ANALYSIS 

o 

o 

<Xv 

0 

0.000372 

0.000371 

-50.2° 

-46.0° 

0.60° 

0.003728 

0.003667 

-48 . 9° 

-46.4° 

1.20° 

0.007214 

0.007091 

-49.9° 

-47 . 6° 

1 . 80° 

0.010212 

0.010050 

-52.1° 

-49.7° 

2.40° 

0.012559 

0.012363 

-55.3° 

-53.0° 

3.00° 

0 . 014134 

0.013921 

-60.2° 

-57.8° 

3.60° 

0 . 014946 

0 . 014720 

-67 .1° 

-64 .7° 

4.20° 

0 . 015168 

0.014929 

-77 .0° 

-74 .5° 

0 

O 

CD 

0.015224 

0.014966 

-90.5° 

-88.1° 

5.40° 

0 . 015831 

0 .015523 

-107.7° 

-105.4° 

6.00° 


0.017739 


0.017334 


-126.1 


-123 . 9 


TABLE  3 


Data  for  6=5  °  ,  Y=4  °  ,  R,=2 ,  R,=400 ,  p, /p,=0 . 6 ,  c,=1500m/s ,  c,=1410m/s 


RECEIVER 

ANGLE 

PRESSURE 

AMPLITUDE 

PHASE 

ANGLE 

IMAGE 

THEORY 

DOUBLET 

ANALYSIS 

IMAGE 

THEORY 

DOUBLET 

ANALYSIS 

0.05° 

0.000087 

0.000083 

56.3° 

78.8° 

0.50° 

0.000818 

0.000817 

62 . 6° 

78 .8° 

1 . 00° 

0.001532 

0.001558 

62.2° 

78 . 9° 

1.50° 

0 . 002151 

0.002158 

63.0° 

79.0° 

2 . 00° 

0.002537 

0.002559 

62.6° 

79.2° 

2.50° 

0.002704 

0 . 002725 

63.3° 

79.3° 

3.00° 

0.002639 

0.002635 

63.0° 

mm 

3.50° 

0 . 002276 

0.002292 

62.4° 

ES9 

4 . 00° 

0.001699 

0 . 001715 

63.0° 

WM 

4.50° 

0.000926 

0.000942 

60.5° 

78 .8° 

5.00° 

0.000025 

0.000032 

-30.6° 

39.8° 

TABLE  4 


Data  for  6=4°  ,  7=1°  ,  R,=5,  R,=12,  D,/D,=0 . 7,  c,=1500m/s,  c,=1420m/s 


RECEIVER 

PRESSURE 

AMPLITUDE 

PHASE 

ANGLE 

ANGLE 

IMAGE 

DOUBLET 

IMAGE 

DOUBLET 

THEORY 

ANALYSIS 

THEORY 

ANALYSIS 

0 

o 

o 

0.040927 

0 . 040984 

-149.3° 

-144.2° 

O 

O 

o 

0.394236 

0.394086 

-148.6° 

-143.5° 

o 

cn 

o 

o 

0.701851 

0.701115 

-146.5° 

-141 .4° 

1.20° 

0 . 873027 

0 . 871301 

-142.3° 

-137 . 3° 

1 . 60° 

0.915180 

0.912006 

-135.2° 

-130.2° 

2 . 00° 

0.887048 

0.882283 

-124 . 6° 

-119.6° 

2 . 40° 

0 . 851987 

0 . 846756 

-112 .8° 

-107.7° 

2 . 80° 

0.810276 

0 . 805379 

-104.5° 

-99.4° 

3.20° 

0 . 704744 

0 . 699685 

-103.6° 

-98 .7° 

3.60° 

0.503485 

0.498563 

-115.2° 

-110 . 9° 

O 

O 

O 

*3* 

0.315934 

0.320342 

-161 .8° 

-158 . 9° 
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TABLE  5 


Data  for  3=3°,  y=1°/Ri=1,  R,=10.  D.,/o,=0 . 9,  c,  =  1500m/s,  c,=1475m/s 


RECEIVER 

PRESSURE 

AMPLITUDE 

ANGLE 

IMAGE 

THEORY 

DOUBLET 

ANALYSIS 

0 . 03° 

0 . 000123 

0 . 000127 

0.30° 

0 . 001255 

0 . 001257 

o 

cr> 

o 

o 

0.002415 

0.002414 

0 . 90° 

0.003352 

0.003377 

1.20° 

0.004028 

0 . 004066 

1.50° 

0 . 004384 

0 . 004421 

0 

O 

00 

iH 

0.004369 

0.004409 

2 . 10° 

0.004010 

0 . 004030 

2.40° 

0.003339 

0 . 003329 

o 

O 

r- 

csi 

0.002471 

0 . 002444 

O 

O 

o 

CO 

0.001898 

0 . 001799 

PHASE  ANGLE 


IMAGE  DOUBLET 


THEORY  ANALYSIS 


-54.4’ 


-49.8° 


-48 .8° 


-49.4’ 


-48.5' 


-47.5° 


-45.7° 


-42 . 1’ 


-35.2° 


-20.8° 


12.7° 


-47 .2° 


-47.2’ 


-47 .0° 


-45.1° 


-43.3° 


-40.2° 


-34 . 1’ 


-20 .1° 


13.5° 


5' 


APPENDIX  "B" 


PROGRAM  URTEXT 

IMAGE  THEORY — CROSS  SLOPE  PRESSURE  VS  RECEIVER  ANGLE (D) 


INTEGER  I , II, N, SI, S2 , N1 , A 

REAL* 4  B,CC,  C2,D1,D2,G,PI,P1,P2,'21,R1,R2,T,T4, 

*  T6, YO, Yl, Y2, Zl, Z2, Z3, Z4, Z5, Z6, T1  (900)  , 

*  R8 (900) , R9 (900) , S (900) ,  C(900) ,E(900) , 

*  F  (900)  , Z, R3, AL, PZ, D 

DOUBLE  PRECISION  Y,W0,W1 

REAL  TQQ, TQQ1 , TQQ2 , TQQ3 

PI  =ACOS (-1 . 0D00) 

WRITE ( * , * ) 'INPUT  B,G,D1,CC,R1,R2,AL,Y0,A' 

READ  (6,  * )  B ,  G ,  D 1 ,  CC ,  R1 ,  R2 ,  AL ,  Y 0 ,  A 
WRITE (6, 2)  B,G,DlrCC,Rl,R2,AL, Y0,A 
2  FORMAT ('B  =  ',F4.1,/,'G  =  ',F4.1,/,'D1  =  ' 

*  ,F4 . 1, /,  ' CC  =  ' ,F4 . 1, /, 'R1  -  ' ,F4 . 1, /, ' R2  = '  , F4 . 1 , /  ,  ' AL 

*  =  ' ,F6.3, */, ' YO  =  ' , F4 . 1 , / , ' A  =  ',14) 

£*******★*************************************************** 
C  INPUT  PARAMETERS 

q*********************************************************** 
C  B  =  WEDGE  ANGLE  (DEG) 

C  G  =  SOURCE  ANGLE  (DEG) 
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C  D  =  RECEIVER  ANGLE  (DEG) 

C  Nl=  #  OF  IMAGE  POINTS  =  INT(180/B) 

C  R1  =  SOURCE  DISTANCE  (IN  DUMP  DISTANCES)  FROM  APEX 

C  R2  =  RECEIVER  DISTANCE  (IN  DUMP  DISTANCES)  FROM  APEX 

C  YO  =  DISTANCE  (IN  DUMP  DISTANCES)  ALONG  APEX 

C  D 1  =  RH01/RH02 

C  CC  =  C1/C2 

C  AL  =  ALPHA/K2 

C  A  =  #  OF  RECEIVER  POSITIONS 

c  - 

C  CHOOSE  SLOW  OR  FAST  BOTTOM  BY  VALUE  OF  SPEED  RATIO 
C  CC  =  C1/C2  . 

C  - 

0********** 

C 

Q********** 

N1  = 

T6  = 

B  = 

G  = 

C2  = 

TQQ 

C - 

C  DECISION  ABOUT  SLOW  OR  FAST  BOTTOM 
c - 


MAIN  PROGRAM 

INT (180. /B) 

180. /PI 
:  B/T6 
B- (G/T6) 

CC*  *2 
=  TAN  (B) 
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C  T4  =  PI/ (2*SIN (ACOS (CC) ) *TAN (B) )  FOR  FAST  BOTTOM 

C  T4  =  PI/ (2*TAN (ACOS (1/CC) ) *TAN (B) )  FOR  SLOW  BOTTOM 

C - 

C 

IF  (CC.LT.l)  THEN 
TQQ1  =ACOS (CC) 

TQQ2  =  SIN(TQQl) 

ELSE 

TQQ1  =  ACOS (1/CC) 

TQQ2  =  TAN(TQQl) 

END  IF 

TQQ3  =  2 . *TQQ2*TQQ 

T4  =  PI/TQQ3 

Q1  =  1/DSQRT (2 . 0D00) 

C 

DO  20  Ml  =  1,  A 
D  =  B- (B/A) *M1 
D2  =  (Y0*Y0) + (R1*R1) + (R2*R2) 

R3  =  2.*R1*R2 
SI  =  1 


C  | THIS  DO  LOOP  CALCULATES  THE  THETA (N)  I 

C  | AND  THE  IMAGE  SLANT  RANGES  R0 (N)  AND  R9 (N) | 


DO  30  N  =  1,  N1 

IF(Sl.GT.O)  T1 (N)  =  (N-1)*B+G 
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IF(Sl.LT.O)  T1 (N)  =  N*B-G 
SI  =  -  SI 

R8 (N)  =  SQRT (D2-R3*COS (Tl (N) -D) ) 
R9(N)  =  SQRT (D2-R3*COS (Tl (N) +D) ) 
30  CONTINUE 


C  | SUM  THE  PRESSURE  OVER  ALL  IMAGES | 

C  | - , 

PI  =  0.0 
P2  =  0.0 

DO  40  N  =  1,  N1 

S2  =  (-1) ** (INT (N/2) ) 


C  | REFLECTION  COEFFICIENTS  ALONG  NTH  UPPER  PATH | 


W1  =  2*C2*AL 
II  =  INT ( (N-l) /2) 

DO  50  I  =  1,11 

S(I)  =  (R1*SIN'T1 (N) -2*I*B) 

&  +R2*SIN (2*I*B-D) ) /R8 (N) 

IF  (S  (I)  .  GE  .  1 )  S(I)=1 
C(I)  =  SQRT (1- (S (I) *S  (I) ) ) 

T  =  S (I) /D1 

W0  =  (-C2+ (C (I) *C  (I) ) ) 

Y  =  DSQRT ( (W0*W0) + (W1*W1) ) 

Z  =  W0 
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IF(Y.LE.Z)  Y  =  Z 
Y1  =  Q1*SQRT (Y+WO) 

Y2  =  — Q1*SQRT (Y— WO) 

Z1  =  T-Y2 
Z2  =  — Y1 

Z3  =  Z1 / (Zl*Zl+Z2*Z2) 

Z4  =  -Z2/ <Z1*Z1+Z2*Z2) 

Z1  =  T+Y2 
Z2  =  Y1 

Z5  =  Z1*Z3— Z2*Z4 

Z6  =  Z1*Z4+Z2*Z3 
E ( I )  =  Z5 
F  ( I )  =  Z6 

50  CONTINUE 

C| - 

C  PRODUCT  OF  REFLECTION  COEFFICIENTS  ALONG  NTH 

C  UPPER  PATH 

C| - 

Z1  =  0 

Z2  =  0 

Z3  =  0 

Z4  =  0 

Z5  =  1 

Z6  =  0 

IF (N. LE . 2 . 00)  GO  TO  110 
DO  60  I  =  1,  II 
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Z1  =  E  ( I ) 


Z2  =  F  ( I ) 

Z3  =  Z5 
Z4  =  Z6 

Z5  =  Z1*Z3-Z2*Z4 
Z6  =  Z1*Z4+Z2*Z3 
60  CONTINUE 

110  Z1  =  Z5 

Z2  =  Z  6 
T  =  T4  *R8 (N) 

Z3  =  COS (T) 

Z4  =  —SIN (T) 

Z5  =  Z1*Z3-Z2*Z4 
Z6  =  Z1*Z4+Z2*Z3 
PI  =  P1+S2*Z5/R8 (N) 
P2  =  P2+S2*Z6/R8 (N) 
11=11+1 


C  | REFLECTION  COEFFICIENTS  ALONG  NTH  LOWER  PATH | 


DO  70  I  =  1,  II 

S(I)  =  (R1*SIN (T1 (N) —2* (I— 1) *B) 

&  +  R2*SIN (2* (1-1) *B+D) ) /R9 (N) 

IF(S(I)  .  GT  .  1 )  S(I)=1 
C (I)  =  SQRT  < 1-S  (I) *S  (I) ) 

T  =  S (I) /D1 
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Cl 

c 

c 

Cl 


WO  =  -C2+C  (I)  *C(I) 

Y  =  DSQRT ( (WO*WO)  +  (W1*W1) ) 
Z  =  WO 

IF(Y.LE.Z)  Y  =  Z 
Y1  =  Q1*SQRT (Y+WO) 

Y2  =  -Q1*SQRT (Y-WO) 

Z1  =  T-Y2 
Z2  =  — Y1 

Z3  =  Zl/  (Z1*Z1+Z2*Z2) 

Z4  =  -  Z2/ <Z1*Z1+Z2*Z2) 

Zl  =  T+Y2 
Z2  =  Y1 

Z5  =  Z1*Z3-Z2*Z4 
Z6  =  Z1*Z4+Z2*Z3 
E ( I )  =  Z5 
F (I)  =  Z6 
CONTINUE 


PRODUCT  OF  REFLECTION  COEFFICIENTS  ALONG  NTH  LOWER 
PATH 


Zl  -  0 
Z2  =  0 
Z3  =  0 
Z4  =  0 
Z5  =  1 
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Z6  =  0 


DO  80 

I 

II 

H 

Z1 

= 

E(I) 

Z2 

= 

F  (I) 

Z3 

= 

Z5 

Z4 

= 

Z6 

Z5 

= 

Z1*Z3— Z2*Z4 

Z6 

= 

Z1*Z4+Z2*Z3 

CONTINUE 

Z1  =  Z5 
Z2  =  Z6 
T  =  T4  *R9 (N) 

Z3  =  COS (T) 

Z4  =  —SIN  (T) 

Z5  =  Z1*Z3-Z2*Z4 
Z6  =  Z1*Z4+Z2*Z3 
PI  =  P1+S2*Z5/R9 (N) 

P2  =  P2+S2*Z6/R9 (N) 

40  CONTINUE 

PZ  =  R1*SQRT <P1*P1+P2*P2) 

PH  =  ATAN2 ( P  2 , P 1 ) *  T  6 
D  =  (B-D) *T6 
WRITE (6, 3)  D,  PZ,  PH 

3  FORMAT (3X, ' D  =  ' , F8 . 4 , 3X, ' PZ  =  ' , F9 . 5, 3X, ' PH  =  ',F9.5) 
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20 


CONTINUE 

END 


66 


APPENDIX  "C" 


PROGRAM  DOUBLT1C 

C*********************************************************** 

C  THIS  PROGRAM  PRODUCES  THE  OUTPUTS  OF  THE  DOUBLET  ANALYSIS 

C  OF  THE  IMAGE  THEORY.  THE  SOURCE  AND  ITS  FIRST  SURFACE 

C  REFLECTED  IMAGE  ARE  TREATED  AS  A  SEPARATE  DOUBLET,  AND 
C  THE  REST  OF  THE  IMAGES  ARE  TREATED  AS  PAIRS  OF  ACOUSTIC 
C  DOUBLETS.  THIS  PROGRAM  IS  INITIALY  MADE  TO  CALCULATE  THE 
C  REFLECTION  COEFFICIENTS  FOR  A  SLOW  BOTTOM  WITHOUT  ANY 
C  ABSORPTION  IN  THE  BOTTOM.  THAT  WAY,  THE  REFLECTION 
C  COEFFICIENTS  ARE  REAL,  AND  THE  PROGRAM  IS  SIMPLER  TO 

C  INITIALY  VERIFY  THE  VALIDITY  OF  THIS  NEW  ANALYSIS.  IT 

C  ALSO  WORKS  FOR  DOWNSLOPE/UPSLOPE  CASE  ONLY.  THE  INITIAL 
C  SETUP  (BEFORE  REACHING  THE  DOLOOP  #20)  IS  APPLICABLE  TO 

C  ALL  OTHER  CASES  (CROSS  SLOPE,  FAST  BOTTOM,  ABSORPTION  JN 

C  THE  BOTTOM) ,  BUT  THE  PROGRAM  WOULD  REQUIRE  MODIFICATIONS 
C  IN  DOLOOP  #20  AND  BEYOND 

INTEGER  A, N1 , NQ1 , Ml , K, M, MQ, N 

DOUBLE  PRECISION  B, G, D1 , Cl , C2 , Y0 , R1 , R2 , R3 , CC, PI , D2 , 
*TQQ, TQQ1,PQ2 (5000) , TQQ2 , TQQ3, T4 , T6, R0P, SNS0P,P1,P2, 

*Z1 (90) ,Z2 (90) , GNP (90) , THP (90) ,EP (90) ,PHIP (90) , 

*PSIP (90) , Z3 (90) , Z4 (90) ,T11 (90) ,  RL (90) , RU (90) , SL, SU, 
*RLD ( 90 ) , RUD (90) ,DRU(90) ,DRL(90) ,SUMRU(90) ,  SUMRL(90) , 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


*Z5 (90) , Z6 (90) , Z7 (90) , Z8 (90) ,RN0P (90) ,D (5000) , PQ1 (5000) 
PI  =  ACOS (-1 . 0D00) 

WRITE (*, *) ' INPUT  B, G, R1,R2, Y0,D1, Cl, C2, A' 

READ (5, *)B,G,R1,R2, Y0 , D1 , Cl , C2 , A 
WRITE ( 6, 2) B, G, Rl , R2  ,Y0,D1,C1,  C2 ,  A 
2  FORMAT ( 9X, ' B  =  ',F5.2,/,  6X,'G  =  ',F6.2,/,  6X,'R1  = 

* ' , F5 .2,1,  6X,'R2  =  ' , F5 . 2 , / ,  6X,'Y0  =  ',F3.1,/, 

*6X,'D1  =  ' , F4 . 2 , / ,  6X, ' Cl  =  ',F7.2,/,  6X,'C2  = 

*' t F10 . 4, /,  6X, / , 9X, ' A  =  ',14) 

'k'k'k’k'k'k'k'k*'k’k'k'k*'k'k'k,k-ki<'k'k'k'k'k'k'k'k’k'k'k'k'k'k'k’k'k'k'k’k'k'k'kit'kitiK*'k'k'k1<1i'k’k**iK*'k* 

* INPUT  PARAMETERS 

************************************************************* 
FR  =  FREQUENCY 
X  =  SCALING  DISTANCE 
B  =  WEDGE  ANGLE (DEG) 

G  =  SOURCE  ANGLE  MEASURED  FROM  THE  SURFACE  (DEG) 

D  =  RECEIVER  ANGLE  MEASURED  FROM  THE  SURFACE  (DEG) 
N1  =  NUMBER  OF  MONOPOLES ( IMAGE  POINTS) 

NQ1  =  NUMBER  OF  PAIRS  OF  ACOUSTIC  DOUBLETS 
Rl  =  SOURCE  DISTANCE (IN  DUMP  DISTANCES)  FROM  APEX 
R2  =  RECEIVER  DISTANCE (IN  DUMP  DISTANCES)  FROM  APEX 
Y0  =  SOURCE -RECEIVER  HORIZONTAL  DISTANCE  IN  DUMP 
DISTANCES  ALONG  APEX 
D1  =  RH01/RH02 

CC  =  C1/C2 

A  =  NUMBER  OF  RECEIVER  POSITIONS 
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C  CHOOSE  SLOW  OR  FAST  BOTTOM  BY  VALUE  OF  SPEED  RATIO 

C  CC  =  C1/C2 

C*********************************************************** 

C*********************************************************** 

C  MAIN  PROGRAM 

C*********************************************************** 

K  =  0 

CC  =  C1/C2 

N1  =  INT (180 . /B) 

NQ1  =  INT ( (INT (360. /B) -2) /4) 

T6  =  180. /PI 
B  =  B/T6 
G  =  G/T6 
D2  =  1./D1 
TQQ  =  TAN (B) 

C 

c  ******************************************************** 

C  DECISION  ABOUT  SLOW  OR  FAST  BOTTOM 

C  T4  =  PI/ (2*SIN (ACOS (CC) ) *TAN (B) )  FOR  FAST  BOTTOM 

C  T4  =  PI/ (2*TAN (ACOS (1/CC) ) *TAN (B) )  FOR  SLOW  BOTTOM 

c  ******************************************************** 
IF(CC.LT.l)  THEN 
TQQ1  =  ACOS (CC) 

TQQ2  =  SIN(TQQl) 

ELSE 
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TQQ1  =  ACOS (1/CC) 

TQQ2  =  TAN(TQQl) 

END  IF 

TQQ3  =  2 . *TQQ2*TQQ 
T4  =  PI/TQQ3 
PRINT*, ' K1X  =  '  ,T4 
C  X  =  T4*C1/ (2*PI*FR) 

CALL  EXCMS ( ' FILEDEF  7  9  CLEAR' ) 

DATA  NOA/79/ 

CALL  EXCMS ('FILEDEF  79  DISK  IFDOUT4  PLOTTER  A1 ' ) 

R3  =  2*R1*R2 
DO  20  Ml  =  1 , A 
K  =  K+l 

C  START  COUNTING  THE  RECEIVERS  POSITIONS 

D (Ml )  =  (B/A) *M1 
DO  N  =  1,N1 

C  CALCULATION  OF  THE  ANGLE  OF  EACH  IMAGE  FROM  THE  BOTTOM 
IF ( (N/2) . LT . (N/2 . ) )  THEN 
Til (N)  =  N*B-G 
ELSE 

Til (N)  =  (N-l) *B+G 
END  IF 

C  CALCULATION  OF  THE  IMAGE  SLANT  RANGES 
C  LOWER  IMAGES 

RLD(N)  =  DSQRT (Rl**2+R2**2fY0**2- 
&  -R3*COS (Til (N) +B-D (Ml) ) ) 
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C  UPPER  IMAGES 

RUD(N)  =  DSQRT (Rl**2+R2**2+Y0**2- 
&  -R3*C0S (Til (N) -B+D (Ml) ) ) 

ENDDO 

C  CALCULATION  OF  THE  REFLECTION  COEFFICIENTS  FOR  A  SLOW 
BOTTOM  WITHOUT  ANY  ABSORPTION 
C  LOWER  IMAGES 

DO  N  =  1 , N1 
RL(N)  =1. 

MQ  =  INT ( (N— 1 ) / 2 ) 

DO  M  =  0,MQ 

SKI  =  R2*DSIN ( (2*M+1) *B-D (Ml) ) 

SL  =  DASIN ( (R1*DSIN (Til (N) -2*M*B) +SK1) /RLD (N) ) 

A50  =  DSQRT (CC**2- (COS (SL) ) **2) 

VL  =  (D2*DSIN(SL) -A50) / (D2*DSIN (SL) +A50) 

RL  (N)  =  RL  (N)  *VL 

ENDDO 
ENDDO 

C  UPPER  IMAGES 

DO  N  =  3 ,  N1 

C  THIS  DOLOOP  STARTS  FROM  THE  THIRD  UPPER  IMAGE  BECAUSE  THE 

C  FIRST  AND  THE  SECOND  UPPER  IMAGES  ARE  CONSIDERED  AS  A 

C  NEUTRAL  DOUBLET 
RU (N)  =  1. 

MQ  =  INT ( (N— 1 ) /2) 

DO  M  =  1 , MQ 
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SK2  =  R2*DSIN ( (2*M~1) *B+D (Ml) ) 

SU  =  DASIN ( (R1*DSIN (Til (N) -2KM*B) +SK2) /RUD (N) ) 
B50  =  DSQRT (CC**2- (COS (SU) ) **2) 

VU  =  (D2*DSIN (SU) -B50) / (D2*DSIN (SU) +B50) 

RU  (N)  =  RU  (N)  *VU 
ENDDO 
ENDDO 

C  SUM  AND  DIFFERENCE  OF  THE  REFLECTION  COEFFICIENTS 

C  BETWEEN  THE  TWO  IMAGES  OF  EACH  ACOUSTIC  DOUBLET 

DO  N  =  1 , NQ1 

IF  ( (2*  (N+2) )  .LE.N1)  THEN 
DRU(N)  =  RU (2*N+1) -RU (2*N+2) 

SUMRU(N)  =  RU  (2*N+1)  +RU  (2*N+2) 

DRL(N)  =  RL(2*N-1) -RL(2*N) 

SUMRL(N)  =  RL (2*N-1) +RL (2*N) 

ELSE 

DRU (N)  =  0 
SUMRU (N)  =  u 
DRL(N)  =  0 
SUMRL(N)  =  0 
END  IF 
ENDDO 

C  PRESSURE  FIELD  FROM  THE  NEUTRAL  DOUBLET  (EQ.31) 

ROP  =  SQRT (R1**2+R2**2+YO**2-2*R1*R2*COS (D (Ml) ) ) 
SNSOP  =  (R2/R0P) *D (Ml) 

PI  =  (2/ROP) *SIN (T4*R1*G*SNS0P) * (SIN (T4*R0P) ) 
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P2  =  (2/ROP) *SIN (T4*R1*G*SNS0P) * (COS (T4*R0P) ) 

C  CALCULATION  OF  THE  TOTAL  PRESSURE  FIELD.  ADD  THE 

C  FIELD  OVER  ALL  DOUBLETS  ACCORDING  TO  EQS . 51  AND  52 

DO  30  N  =  1 , NQ1 
C  EQ.35 

RNOP (N)  =  SQRT (Rl**2+R2**2+Y0**2- 
&  — 2*Rl*R2*COS (2*N*B) ) 

C  EQ.45 

GNP (N)  =  T4*R1*R2/RN0P (N) 

C  EQ.46 

THP (N)  =  GNP  ;N) *G*SIN (2*N*B) 

C  EQ.48 

EP (N)  =  GNP (N) *G*D (Ml ) *COS (2*N*B) 

C  EQ.47 

PHIP (N)  =  GNP <N) *D<M1) *SIN(2*N*B) 

C  EQ.49 

KSIP(N)  =  0 . 3* (D (Ml) *EP (N) /G— 

&  -PHIP <N) **2/ (T4*RN0P (N) ) ) 

C  APPLY  EQ  32 

Z1  <N)  =  SIN  (THP  (N)  +EP  (N)  )  * 

&  *SIN (T4*RN0P (N) +PHIP (N) +KSIP (N) ) 

Z1 (N)  =  SUMRU (N) *Z1 (N) 

Z2 (N)  =  COS (THP (N) +EP (N) ) * 

&  *COS (T4*RN0P (N) +PHIP (N) +KSIP (N) ) 

Z2 (N)  =  DRU (N) *Z2 (N) 

Z3(N)  =  SIN (THP (N) -EP (N) ) * 
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&  *  SIN (T4*RN0P (N) -PHIP (N) +KSIP (N)  ) 

Z3(N)  =  -SUMRL (N) *Z3 (N) 

Z4  (N)  =  COS (THP (N) -EP  (N)  ) * 

&  *COS (T4*RN0P (N) -PHIP (N) +KSIP  (N) ) 

Z4 (N)  =  -DRL (N) *Z4 (N) 

Z5(N)  =  SIN  (THP  (N)  +EP  (N)  )  * 

&  *COS (T4*RN0P (N) +PHIP (N) +KSIP (N) ) 

Z5{N)  =  SUMRU(N)  *Z5  (N) 

Z6(N)  =  COS (THP (N) +EP (N) ) * 

&  *SIN (T4*RN0P (N) +PHIP (N) +KSIP (N) ) 

Z6(N)  =  -DRU(N)  *Z6(N) 

Z7  (N)  =  COS  (THP  (N) -EP  (N)  )  * 

&  *SIN (T4*RN0P (N) -PHIP (N) +K .  "  (W) ) 

Z7 (N)  =  DRL (N) *Z7 (N) 

Z8  (N)  =  SIN  (THP  (N)  -EP  (N)  ) 

&  *COS (T4*RN0P (N) -PHIP (N) +KSIP (N) ) 

Z8 (N)  =  -SUMRL (N) *Z8 (N) 

PI  =  (-1) **N* (1 . /RNOP (N) ) * 

&  * (Z1 (N) +Z2 (N) +Z3 (N) +Z4 (N) ) +P1 

P2  =  (-1) **N* (1 . /RNOP  (N) ) * 

&  * (Z5 (N) +Z6 (N) +Z7 (N) +Z8 (N) ) +P2 

30  CONTINUE 
C  AMPLITUDE 

PQl(Ml)  =  R1*SQRT(P1**2+P2**2) 

C  PHASE 

PQ2 (Ml )  =  ATAN2 (P2 , PI ) *T6 
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C  CONVERT  THE  RECEIVER  ANGLE  IN  DEGREES 

D (Ml )  =  D (Ml) *180 . /PI 
C  WRITE (6, 3)D (Ml) , PQ1 (Ml) ,  PQ2 (Ml) 

C  3  FORMAT (9X, F6 . 3, 8X, F9.6,8X,F10.5) 

20  CONTINUE 

C  STORE  THE  RESULTS  IN  THE  FILE  IFDOUT4  PLOTTER  A1 .  FROM 

C  THAT  FILE  THEY  ARE  RECALLED  IN  URTEXT1C  EXECUTION,  FOR 

C  COMPARISON  WITH  THE  IMAGE  THEORY  RESULTS . 

DO  21  J  =  1 , K 

WRITE (NOA, * ) D  ( J)  , PQ1 (J) , PQ2 (J) 

21  CONTINUE 
END 
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